% L1 and L2 metric of P
clc
clear
load ('duration1yrs_TP06YRJfinR5G2501000kMR7.mat','durationCount', 'durationProb', 'NUMBER_OF_YEAR','DEFAULT');



RATE_LIST= 1.5:0.5:4.5;
RATE_LIST=[RATE_LIST, DEFAULT];
RATE_LIST_LENGTH = length(RATE_LIST);
X = (1:1:12);


avgProb = zeros(RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
L1 = zeros(NUMBER_OF_YEAR,1);
L2 = zeros(NUMBER_OF_YEAR,1);
WAD = zeros(NUMBER_OF_YEAR,1);
NAD = zeros(NUMBER_OF_YEAR,1);
WSD = zeros(NUMBER_OF_YEAR,1);
NSD = zeros(NUMBER_OF_YEAR,1);
SVD = zeros(NUMBER_OF_YEAR,1);
D1 = zeros(NUMBER_OF_YEAR,1);
D2 = zeros(NUMBER_OF_YEAR,1);
D3 = zeros(NUMBER_OF_YEAR,1);
D4 = zeros(NUMBER_OF_YEAR,1);
D1ns = zeros(NUMBER_OF_YEAR,1);
D2ns = zeros(NUMBER_OF_YEAR,1);
D3ns = zeros(NUMBER_OF_YEAR,1);
D4ns = zeros(NUMBER_OF_YEAR,1);


%  aggregate average transition probability
ss(:, :) = sum(durationCount,1);
ttlCol = sum(ss, 2);
for i=1:RATE_LIST_LENGTH
    avgProb(:, i)= ss(:, i)./ttlCol;
end

% calculate metric distance L1 , L2 , WAD , SVD
for i = 1: NUMBER_OF_YEAR
    p(:, :) = durationProb (i, :, :);
    L1(i) = sum(sum(abs(p - avgProb))); 
    L2(i) = sqrt(sum(sum((p - avgProb).^2)));
    WAD(i) = sum(sum(abs(p - avgProb).*p));
    WSD(i) = sqrt(sum(sum((p - avgProb).^2.*p)));
    SVD(i) = mean(svd(p))- mean(svd(avgProb)); % maybe wrong! - Mar13,2012
end

% calculate metric distance NAD, D2
for i = 1 : NUMBER_OF_YEAR
    p(:, :) = durationProb (i, :, :);
    for j = 1 : RATE_LIST_LENGTH-1
        for k = 1 : RATE_LIST_LENGTH
            
            if p(j , k) ~=0
                mNAD = abs(p(j, k) - avgProb(j, k))/p(j, k);
                mNSD = (p(j, k) - avgProb(j, k))^2/p(j, k);
                mD2 = (j - k)*(p(j, k) - avgProb(j, k))/p(j, k);
                mD4 = (j - k)* sign( p(j, k) - avgProb(j, k) ) * (p(j, k) - avgProb(j, k))^2 / p(j, k);
                NAD(i) = NAD (i)+ mNAD;
                NSD(i) = NSD (i)+ mNSD;           
            end
            
            mD1 = (j - k)*(p(j, k) - avgProb(j, k));
            mD3 = (j - k)* sign(p(j, k) - avgProb(j, k))*(p(j, k) - avgProb(j, k))^2;
            
            if k~= RATE_LIST_LENGTH
                D1(i) = D1(i) + mD1;
                D2(i) = D2(i) + mD2;
                D3(i) = D3(i) + mD3;
                D4(i) = D4(i) + mD4;
                D1ns(i) = D1(i) + mD1;
            else
                D1(i) = D1(i) + mD1*RATE_LIST_LENGTH;
                D2(i) = D2(i) + mD2*RATE_LIST_LENGTH;
                D3(i) = D3(i) + mD3*RATE_LIST_LENGTH;
                D4(i) = D4(i) + mD4*RATE_LIST_LENGTH;
                D1ns(i) = D1(i) + mD1*(RATE_LIST_LENGTH)^2;
                D2ns(i) = D2(i) + mD2*(RATE_LIST_LENGTH)^2;
                D3ns(i) = D3(i) + mD3*(RATE_LIST_LENGTH)^2;
                D4ns(i) = D4(i) + mD4*(RATE_LIST_LENGTH)^2;
            end
            
        end
    end
end

NSD = sqrt(NSD);

figure;
subplot(1,2,1); plot (X,L1,'-*',X,L2,':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('L1, L2 Norms');
legend('L1','L2','Location','NorthWest');
subplot(1,2,2); plot (X,zscore(L1),'-*', X,zscore(L2),':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized L1, L2 Norms');
legend('L1','L2','Location','NorthWest');

figure;
subplot(1,2,1); plot (X,WAD,'-*', X,NAD,':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('WAD, NAD ');
legend('WAD','NAD','Location','NorthWest');
subplot(1,2,2); plot (X,zscore(WAD),'-*', X,zscore(NAD),':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized WAD, NAD ');
legend('WAD','NAD','Location','NorthWest');

figure;
subplot(1,2,1); plot (X,WSD,'-*', X,NSD,':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('WSD, NSD ');
legend('WSD','NSD','Location','NorthWest');
subplot(1,2,2); plot (X,zscore(WSD),'-*', X,zscore(NSD),':o')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized WSD, NSD ');
legend('WSD','NSD','Location','NorthWest');

figure;
subplot(1,2,1);plot (X,SVD,'-*')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('SVD ');
legend('SVD','Location','NorthWest');
subplot(1,2,2); plot (X,zscore(SVD),'-*');
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized SVD ');
legend('SVD','Location','NorthWest');

figure;
subplot(1,2,1); plot (X,D1,'-*', X,D2,'-o', X,D3,'-+',X,D4,'-s')
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('D1, D2, D3, D4');
legend('D1', 'D2', 'D3','D4','Location','NorthWest');
subplot(1,2,2); plot (X,zscore(D1),'-*', X,zscore(D2),'-o', X,zscore(D3),'-+',X,zscore(D4),'-s');
set(gca,'XTick',1:1:12);
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized D1, D2, D3 D4');
legend('D1', 'D2', 'D3','D4','Location','NorthWest');

figure;
subplot(4,2,1); plot (X,D1,'-*', X,D1ns,'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('D1, D1ns');
legend('D1', 'D1ns','Location','NorthWest');
subplot(4,2,2); plot (X,zscore(D1),'-*', X,zscore(D1ns),'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized D1, D1ns');
legend('D1', 'D1ns','Location','NorthWest');
subplot(4,2,3); plot (X,D2,'-*', X,D2ns,'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('D2, D2ns');
legend('D2', 'D2ns','Location','NorthWest');
subplot(4,2,4); plot (X,zscore(D2),'-*', X,zscore(D1ns),'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized D2, D2ns');
legend('D2', 'D2ns','Location','NorthWest');
subplot(4,2,5); plot (X,D3,'-*', X,D3ns,'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('D3, D3ns');
legend('D3', 'D3ns','Location','NorthWest');
subplot(4,2,6); plot (X,zscore(D3),'-*', X,zscore(D3ns),'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized D3, D3ns');
legend('D3', 'D3ns','Location','NorthWest');
subplot(4,2,7); plot (X,D4,'-*', X,D4ns,'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('D4, D4ns');
legend('D4', 'D4ns','Location','NorthWest');
subplot(4,2,8); plot (X,zscore(D4),'-*', X,zscore(D4ns),'--*')
set(gca,'XTick',1:1:12)
set(gca,'XTickLabel',{'1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'});
xlabel('Year');
ylabel('Standadized D4, D4ns');
legend('D4', 'D4ns','Location','NorthWest');

